jbagnall@broadinstitute.org, adapted from scripts by Marek Orzechowski
Jan. 21, 2025
Calculate moderated z-score based on RNAseq data. This notebook is to be run after 01_qc_and_vst.Rmd. This notebook performs the points below:
  1. Quantile normalize
  2. Robust z-score
  3. Moderated z-score, collapsing replicates based on correlation between replicates

User inputs

#working directory
wkdir = '/Users/sophie_chen/downloads/psa_rnaseq_manuscript_rproject-main/'

#output directory
outdir = '/Users/sophie_chen/downloads/psa_rnaseq_manuscript_rproject-main/output'

#run_name is something to identify this data set, will be included in saved file names
run_name = "mocp_0193_6h"

#run_name_suffix is something to identify this data set (in case it's a subset of run_name), will be included in saved file names, can be empty string
run_name_suffix = ""

#project_id0 could be the same as run_name, but is something corresponding to a field in the input data frame to identify this data set
project_id0 = "mocp_0193"

#date
date0 = "25127"

#vst file path, VST data calculated in 01_qc_and_vst.Rmd script
#vst_filepath = paste0(outdir, paste(run_name, run_name_suffix, "Project_vst_remove_low_count_samples_and_genes_n120x5849.gct", sep = "_"))
vst_filepath = "/Users/sophie_chen/downloads/psa_rnaseq_manuscript_rproject-main/output/mocp_0193_6h_vst_remove_low_count_samples_and_genes_n180x5846.gct"

#Column name in the column metadata from the VST file that identifies the batches within which you will calculate z-scores, e.g. something that identifies the 384 well plates, or different data sets
batch_column_name = "project_plate_id_384well"

#Value at which to clip the absolute value of the resulting z-scores, so outliers don't throw off down-stream correlations to the reference set
#This value is usually 10 for z-scores (larger # of samples), and larger, e.g. 500, for nz-scores (smaller # of samples)
clip_value = 10
                     

Load libraries & data paths

Quantile normalization

Group data by normalization units (comb_id)

library(cmapR)
library(dplyr)
source('./functions/Functions_general.R')
source('./functions/Functions_mod_z_score.R')

end_char = substr(outdir, nchar(outdir), nchar(outdir))
if(end_char != "/"){
  outdir = paste0(outdir, '/')
}

if(!dir.exists(outdir)){
  print("Creating output directory")
  dir.create(outdir)
}

getwd()
[1] "/Users/sophie_chen/Downloads/psa_rnaseq_manuscript_rproject-main"
#VST data calculated in 01_qc_and_vst.Rmd script
data_vst = parse_gctx(vst_filepath)
parsing as GCT v1.3
/Users/sophie_chen/downloads/psa_rnaseq_manuscript_rproject-main/output/mocp_0193_6h_vst_remove_low_count_samples_and_genes_n180x5846.gct 5846 rows, 180 cols, 16 row descriptors, 25 col descriptors
data_mat = data_vst@mat
row_meta = data_vst@rdesc
col_meta = data_vst@cdesc

#Compound metadata
# comp_meta = read.csv('./reference_files/Metadata_compound_231129.csv', stringsAsFactors = F)

Quantile normalization within each normalization unit (normalization across each sample)

#Define groups (batches) to do quantile normalization
temp = col_meta[[batch_column_name]]
col_meta[["comb_id"]] = temp
remove(temp)

comb_id_unique = unique(col_meta$comb_id)
print(comb_id_unique)
[1] "MOCP_0193_1"
#update metadata in gctx file
data_vst@cdesc = col_meta

Plotting qnorm data for QC visualization

#Loop over each comb_id

data_qnorm = data_mat
flag0 = rep(0, length(comb_id_unique)) #flag to check if all comb_ids ran
for(num in 1:length(comb_id_unique)){
    idx = which(col_meta$comb_id == comb_id_unique[num])
    sum(length(idx))
    if(sum(length(idx))>0){
      data_qnorm[,idx] = qnorm_median(data_mat[,idx])
    }else{
      flag0[num] = 1
    }
}

#If any comb_ids had no samples, sum(flag0) > 0
if(sum(flag0) > 0){
  print("Warning: Not all comb_ids were quantile normalized")
}else{
  print("All comb_ids quantile normalized")
}
[1] "All comb_ids quantile normalized"
any(is.na(data_qnorm))
[1] FALSE
#Make boxplot of each sample
col_meta_select = dplyr::select(col_meta, id, comb_id) %>% distinct()
data_qnorm_df = as.data.frame(data_qnorm)
data_qnorm_df$gene_id = rownames(data_qnorm_df)
data_qnorm_df = tidyr::gather(data_qnorm_df, key = "id", value = "qnorm_vst", -gene_id)
data_qnorm_df = dplyr::left_join(data_qnorm_df, col_meta_select, by = "id")

ggplot(data_qnorm_df, aes(x = id, y = qnorm_vst, group = id)) +
  geom_boxplot()+
  facet_wrap(~comb_id, scales = "free_x")+
  theme_bw(base_size = 14)


data_qnorm_pca = run_pca(data_qnorm, col_meta = col_meta, color_col = "strain_id", scale0 = T)

Z-score across each gene

Prepare with QC

#Save qnorm data to gct
qnorm_file_path = paste0(outdir, paste(run_name,run_name_suffix, "qnorm_remove_low_count_samples_and_genes", sep = "_"))

save_to_gct(data_mat = data_qnorm, data_cdesc = col_meta, data_rdesc = row_meta, savefilepath = qnorm_file_path)
Saving file to  /Users/sophie_chen/downloads/psa_rnaseq_manuscript_rproject-main/output/mocp_0193_6h__qnorm_remove_low_count_samples_and_genes_n180x5846.gct 
Dimensions of matrix: [5846x180]
Setting precision to 4
Saved.

Perform Z-score calculation per row

#Each row is a gene, visualize metrics that will be used for z-score calculation
row_medians = rowMedians(data_qnorm, na.rm = T)
plot(x = seq(1:length(row_medians)), y = row_medians)


row_mads = apply(data_qnorm, MARGIN = 1, FUN = function(x){mad(x, na.rm = T)})
plot(x = seq(1:length(row_mads)), y = row_mads)

Clip extreme z-score values

zscore_list = list()
data_zscore = data_qnorm

flag0 = rep(0, length(comb_id_unique)) #flag to check if all comb_ids ran
for(num in 1:length(comb_id_unique)){
    idx = which(col_meta$comb_id == comb_id_unique[num])
    sum(length(idx))
    if(sum(length(idx))>0){
       zscore_list[[num]] = robust_zscore(data_qnorm[,idx], dim0 = 1) #z-score along rows, outputs zscores, medians and mads
       data_zscore[,idx] = zscore_list[[num]][[1]] #only taking the z-scores
    }else{
      flag0[num] = 1
    }
}

#Any comb_ids didn't get z-scored?
if(sum(flag0) > 0){
  print("Warning: Not all comb_ids were standardized")
}

Plot PCA

#Clip z-score values to account for extreme cases where genes with low count had low variation (i.e. low median absolute deviation)
#Look at genes' z-score distribution
hist(data_zscore, breaks = 100)

quantile(data_zscore, c(0.001, 0.005, 0.01, 0.05, 0.95, 0.99, 0.995, 0.999))
     0.1%      0.5%        1%        5%       95%       99%     99.5%     99.9% 
-5.069348 -3.540297 -2.992669 -1.863136  1.785939  3.143715  3.930308  6.602259 
data_zscore = base::pmax(data_zscore, -1*clip_value)
data_zscore = base::pmin(data_zscore, clip_value)

Save z-scores into GCT file

data_zscore_pca = run_pca(data_zscore, col_meta = col_meta, color_col = "strain_id", scale0 = T)

Collapse replicates

if(!all(colnames(data_zscore) == col_meta$id)){
  print("Warning: column names of matrix & metadata IDs do not match")
}
if(!all(rownames(data_zscore) == row_meta$gene_id)){
  print("Warning: row names do not match gene IDs in metadata")
}

save_to_gct(data_mat = data_zscore, data_cdesc = col_meta, data_rdesc = row_meta, savefilepath = paste0(outdir, paste(run_name,run_name_suffix, "zscore_remove_low_count_samples_and_genes", sep = "_")))
Saving file to  /Users/sophie_chen/downloads/psa_rnaseq_manuscript_rproject-main/output/mocp_0193_6h__zscore_remove_low_count_samples_and_genes_n180x5846.gct 
Dimensions of matrix: [5846x180]
Setting precision to 4
Saved.
num_row_zscore = dim(data_zscore)[1]
num_col_zscore = dim(data_zscore)[2]

Save replicate collapsed data as gct file

zscore = cmapR::parse_gctx(paste0(outdir, paste(run_name,run_name_suffix, "zscore_remove_low_count_samples_and_genes", sep = "_"), "_n", num_col_zscore, "x", num_row_zscore, ".gct"))
parsing as GCT v1.3
/Users/sophie_chen/downloads/psa_rnaseq_manuscript_rproject-main/output/mocp_0193_6h__zscore_remove_low_count_samples_and_genes_n180x5846.gct 5846 rows, 180 cols, 16 row descriptors, 26 col descriptors
col_meta = zscore@cdesc
data_zscore = zscore@mat


#identify replicates
#perform moderated collapsing of averages
condition_id_unique = unique(col_meta$condition_id)
num_genes = dim(data_zscore)[1]
num_samples = dim(data_zscore)[2]
num_conditions = length(condition_id_unique)

data_rep_collapse = matrix(NA, nrow = num_genes, ncol = num_conditions)

all(col_meta$id == colnames(data_zscore))
[1] TRUE
for(num in 1:num_conditions){
  idx = which(col_meta$condition_id == condition_id_unique[num])
  cidx = which(colnames(data_zscore) %in% col_meta$id[idx])
  data_rep_collapse[,num] = modzs_collapse_avg(data = data_zscore[,cidx], ridx = seq(1, num_genes, by = 1))
}
colnames(data_rep_collapse) = condition_id_unique
rownames(data_rep_collapse) = rownames(data_zscore)

Plots and QC

#Revamp column metadata since replicates are collapsed
col_meta_collapse = select(sample_meta, condition_id, pert_id, pert_dose, pert_idose, pert_dose_unit, pert_iname, pert_itime, pert_time, pert_time_unit, pert_type, project_id, strain_id, x_dose_mic_high_inoculum) %>% distinct()

col_meta_collapse_sub = data.frame(condition_id = condition_id_unique)
col_meta_collapse_sub = left_join(col_meta_collapse_sub, col_meta_collapse, by = "condition_id")
all(colnames(data_rep_collapse) == col_meta_collapse_sub$condition_id)
Warning: longer object length is not a multiple of shorter object length
[1] FALSE
all(rownames(data_rep_collapse) == row_meta$gene_id)
[1] TRUE
save_to_gct(data_mat = data_rep_collapse, data_cdesc = col_meta_collapse_sub, data_rdesc = row_meta, savefilepath = paste0(outdir, paste(run_name,run_name_suffix, "modzscore_remove_low_count_samples_and_genes_rep_collapse", sep = "_")))
Error in methods::validObject(ds) : 
  invalid class “GCT” object: cdesc must either have 0 rows or the same number of rows as matrix has columns

Add inactive info and dose ranks to the column metadata if desired (optional)

col_meta_collapse_sub$id = col_meta_collapse_sub$condition_id
data_zscore_rep_collapse_pca = run_pca(data_rep_collapse, col_meta = col_meta_collapse_sub, color_col = "pert_dose", scale0 = T)

LS0tCnRpdGxlOiAiQ2FsY3VsYXRlIG1vZGVyYXRlZCB6LXNjb3JlIGFjcm9zcyBhbGwgdHJlYXRtZW50cywgYnkgYmF0Y2giCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCgpqYmFnbmFsbEBicm9hZGluc3RpdHV0ZS5vcmcsIGFkYXB0ZWQgZnJvbSBzY3JpcHRzIGJ5IE1hcmVrIE9yemVjaG93c2tpPC9icj4KSmFuLiAyMSwgMjAyNTwvYnI+CkNhbGN1bGF0ZSBtb2RlcmF0ZWQgei1zY29yZSBiYXNlZCBvbiBSTkFzZXEgZGF0YS4gVGhpcyBub3RlYm9vayBpcyB0byBiZSBydW4gYWZ0ZXIgMDFfcWNfYW5kX3ZzdC5SbWQuIFRoaXMgbm90ZWJvb2sgcGVyZm9ybXMgdGhlIHBvaW50cyBiZWxvdzogPGJyIC8+CjxvbD4KICA8bGk+UXVhbnRpbGUgbm9ybWFsaXplPC9saT4KICA8bGk+Um9idXN0IHotc2NvcmU8L2xpPgogIDxsaT5Nb2RlcmF0ZWQgei1zY29yZSwgY29sbGFwc2luZyByZXBsaWNhdGVzIGJhc2VkIG9uIGNvcnJlbGF0aW9uIGJldHdlZW4gcmVwbGljYXRlczwvbGk+PC9vbD4KCgojIFVzZXIgaW5wdXRzCmBgYHtyfQojd29ya2luZyBkaXJlY3RvcnkKd2tkaXIgPSAnL1VzZXJzL3NvcGhpZV9jaGVuL2Rvd25sb2Fkcy9wc2Ffcm5hc2VxX21hbnVzY3JpcHRfcnByb2plY3QtbWFpbi8nCgojb3V0cHV0IGRpcmVjdG9yeQpvdXRkaXIgPSAnL1VzZXJzL3NvcGhpZV9jaGVuL2Rvd25sb2Fkcy9wc2Ffcm5hc2VxX21hbnVzY3JpcHRfcnByb2plY3QtbWFpbi9vdXRwdXQnCgojcnVuX25hbWUgaXMgc29tZXRoaW5nIHRvIGlkZW50aWZ5IHRoaXMgZGF0YSBzZXQsIHdpbGwgYmUgaW5jbHVkZWQgaW4gc2F2ZWQgZmlsZSBuYW1lcwpydW5fbmFtZSA9ICJtb2NwXzAxOTNfNmgiCgojcnVuX25hbWVfc3VmZml4IGlzIHNvbWV0aGluZyB0byBpZGVudGlmeSB0aGlzIGRhdGEgc2V0IChpbiBjYXNlIGl0J3MgYSBzdWJzZXQgb2YgcnVuX25hbWUpLCB3aWxsIGJlIGluY2x1ZGVkIGluIHNhdmVkIGZpbGUgbmFtZXMsIGNhbiBiZSBlbXB0eSBzdHJpbmcKcnVuX25hbWVfc3VmZml4ID0gIiIKCiNwcm9qZWN0X2lkMCBjb3VsZCBiZSB0aGUgc2FtZSBhcyBydW5fbmFtZSwgYnV0IGlzIHNvbWV0aGluZyBjb3JyZXNwb25kaW5nIHRvIGEgZmllbGQgaW4gdGhlIGlucHV0IGRhdGEgZnJhbWUgdG8gaWRlbnRpZnkgdGhpcyBkYXRhIHNldApwcm9qZWN0X2lkMCA9ICJtb2NwXzAxOTMiCgojZGF0ZQpkYXRlMCA9ICIyNTEyNyIKCiN2c3QgZmlsZSBwYXRoLCBWU1QgZGF0YSBjYWxjdWxhdGVkIGluIDAxX3FjX2FuZF92c3QuUm1kIHNjcmlwdAojdnN0X2ZpbGVwYXRoID0gcGFzdGUwKG91dGRpciwgcGFzdGUocnVuX25hbWUsIHJ1bl9uYW1lX3N1ZmZpeCwgIlByb2plY3RfdnN0X3JlbW92ZV9sb3dfY291bnRfc2FtcGxlc19hbmRfZ2VuZXNfbjEyMHg1ODQ5LmdjdCIsIHNlcCA9ICJfIikpCnZzdF9maWxlcGF0aCA9ICIvVXNlcnMvc29waGllX2NoZW4vZG93bmxvYWRzL3BzYV9ybmFzZXFfbWFudXNjcmlwdF9ycHJvamVjdC1tYWluL291dHB1dC9tb2NwXzAxOTNfNmhfdnN0X3JlbW92ZV9sb3dfY291bnRfc2FtcGxlc19hbmRfZ2VuZXNfbjE4MHg1ODQ2LmdjdCIKCiNDb2x1bW4gbmFtZSBpbiB0aGUgY29sdW1uIG1ldGFkYXRhIGZyb20gdGhlIFZTVCBmaWxlIHRoYXQgaWRlbnRpZmllcyB0aGUgYmF0Y2hlcyB3aXRoaW4gd2hpY2ggeW91IHdpbGwgY2FsY3VsYXRlIHotc2NvcmVzLCBlLmcuIHNvbWV0aGluZyB0aGF0IGlkZW50aWZpZXMgdGhlIDM4NCB3ZWxsIHBsYXRlcywgb3IgZGlmZmVyZW50IGRhdGEgc2V0cwpiYXRjaF9jb2x1bW5fbmFtZSA9ICJwcm9qZWN0X3BsYXRlX2lkXzM4NHdlbGwiCgojVmFsdWUgYXQgd2hpY2ggdG8gY2xpcCB0aGUgYWJzb2x1dGUgdmFsdWUgb2YgdGhlIHJlc3VsdGluZyB6LXNjb3Jlcywgc28gb3V0bGllcnMgZG9uJ3QgdGhyb3cgb2ZmIGRvd24tc3RyZWFtIGNvcnJlbGF0aW9ucyB0byB0aGUgcmVmZXJlbmNlIHNldAojVGhpcyB2YWx1ZSBpcyB1c3VhbGx5IDEwIGZvciB6LXNjb3JlcyAobGFyZ2VyICMgb2Ygc2FtcGxlcyksIGFuZCBsYXJnZXIsIGUuZy4gNTAwLCBmb3Igbnotc2NvcmVzIChzbWFsbGVyICMgb2Ygc2FtcGxlcykKY2xpcF92YWx1ZSA9IDEwCiAgICAgICAgICAgICAgICAgICAgIApgYGAKCgpgYGB7ciwgc2V0dXAsIGluY2x1ZGU9RkFMU0V9CiMgTWFrZSBzdXJlIGRpcmVjdG9yaWVzIGV4aXN0LCBvdGhlcndpc2UgY3JlYXRlIHRoZW0KZW5kX2NoYXIgPSBzdWJzdHIod2tkaXIsIG5jaGFyKHdrZGlyKSwgbmNoYXIod2tkaXIpKQppZihlbmRfY2hhciAhPSAiLyIpewogIHdrZGlyID0gcGFzdGUwKHdrZGlyLCAnLycpCn0KCmlmKGRpci5leGlzdHMod2tkaXIpKXsKICBzZXR3ZCh3a2RpcikKICBrbml0cjo6b3B0c19rbml0JHNldChyb290LmRpciA9IHdrZGlyKSAjVGhpcyBuZWVkcyB0byBiZSBpbiB0aGUgc2V0dXAgY2h1bmsgdG8gbWFrZSBpdCBnbG9iYWwgZm9yIGFsbCBjaHVua3MKfWVsc2V7CiAgc3RvcCgiV29ya2luZyBkaXJlY3RvcnkgZG9lcyBub3QgZXhpc3QiKQp9CnByaW50KHBhc3RlMCgiV29ya2luZyBkaXJlY3Rvcnk6ICIsIGdldHdkKCkpKQpgYGAKCgoKIyMgTG9hZCBsaWJyYXJpZXMgJiBkYXRhIHBhdGhzCmBgYHtyLCBpbmNsdWRlPUZBTFNFfQpsaWJyYXJ5KGNtYXBSKQpsaWJyYXJ5KGRwbHlyKQpzb3VyY2UoJy4vZnVuY3Rpb25zL0Z1bmN0aW9uc19nZW5lcmFsLlInKQpzb3VyY2UoJy4vZnVuY3Rpb25zL0Z1bmN0aW9uc19tb2Rfel9zY29yZS5SJykKCmVuZF9jaGFyID0gc3Vic3RyKG91dGRpciwgbmNoYXIob3V0ZGlyKSwgbmNoYXIob3V0ZGlyKSkKaWYoZW5kX2NoYXIgIT0gIi8iKXsKICBvdXRkaXIgPSBwYXN0ZTAob3V0ZGlyLCAnLycpCn0KCmlmKCFkaXIuZXhpc3RzKG91dGRpcikpewogIHByaW50KCJDcmVhdGluZyBvdXRwdXQgZGlyZWN0b3J5IikKICBkaXIuY3JlYXRlKG91dGRpcikKfQoKZ2V0d2QoKQoKI1ZTVCBkYXRhIGNhbGN1bGF0ZWQgaW4gMDFfcWNfYW5kX3ZzdC5SbWQgc2NyaXB0CmRhdGFfdnN0ID0gcGFyc2VfZ2N0eCh2c3RfZmlsZXBhdGgpCmRhdGFfbWF0ID0gZGF0YV92c3RAbWF0CnJvd19tZXRhID0gZGF0YV92c3RAcmRlc2MKY29sX21ldGEgPSBkYXRhX3ZzdEBjZGVzYwoKI0NvbXBvdW5kIG1ldGFkYXRhCiMgY29tcF9tZXRhID0gcmVhZC5jc3YoJy4vcmVmZXJlbmNlX2ZpbGVzL01ldGFkYXRhX2NvbXBvdW5kXzIzMTEyOS5jc3YnLCBzdHJpbmdzQXNGYWN0b3JzID0gRikKCmBgYAoKCiMjIFF1YW50aWxlIG5vcm1hbGl6YXRpb24KIyMjIEdyb3VwIGRhdGEgYnkgbm9ybWFsaXphdGlvbiB1bml0cyAoY29tYl9pZCkKYGBge3J9CiNEZWZpbmUgZ3JvdXBzIChiYXRjaGVzKSB0byBkbyBxdWFudGlsZSBub3JtYWxpemF0aW9uCnRlbXAgPSBjb2xfbWV0YVtbYmF0Y2hfY29sdW1uX25hbWVdXQpjb2xfbWV0YVtbImNvbWJfaWQiXV0gPSB0ZW1wCnJlbW92ZSh0ZW1wKQoKY29tYl9pZF91bmlxdWUgPSB1bmlxdWUoY29sX21ldGEkY29tYl9pZCkKcHJpbnQoY29tYl9pZF91bmlxdWUpCgojdXBkYXRlIG1ldGFkYXRhIGluIGdjdHggZmlsZQpkYXRhX3ZzdEBjZGVzYyA9IGNvbF9tZXRhCmBgYAoKIyMjIFF1YW50aWxlIG5vcm1hbGl6YXRpb24gd2l0aGluIGVhY2ggbm9ybWFsaXphdGlvbiB1bml0IChub3JtYWxpemF0aW9uIGFjcm9zcyBlYWNoIHNhbXBsZSkKYGBge3J9CiNMb29wIG92ZXIgZWFjaCBjb21iX2lkCgpkYXRhX3Fub3JtID0gZGF0YV9tYXQKZmxhZzAgPSByZXAoMCwgbGVuZ3RoKGNvbWJfaWRfdW5pcXVlKSkgI2ZsYWcgdG8gY2hlY2sgaWYgYWxsIGNvbWJfaWRzIHJhbgpmb3IobnVtIGluIDE6bGVuZ3RoKGNvbWJfaWRfdW5pcXVlKSl7CiAgICBpZHggPSB3aGljaChjb2xfbWV0YSRjb21iX2lkID09IGNvbWJfaWRfdW5pcXVlW251bV0pCiAgICBzdW0obGVuZ3RoKGlkeCkpCiAgICBpZihzdW0obGVuZ3RoKGlkeCkpPjApewogICAgICBkYXRhX3Fub3JtWyxpZHhdID0gcW5vcm1fbWVkaWFuKGRhdGFfbWF0WyxpZHhdKQogICAgfWVsc2V7CiAgICAgIGZsYWcwW251bV0gPSAxCiAgICB9Cn0KCiNJZiBhbnkgY29tYl9pZHMgaGFkIG5vIHNhbXBsZXMsIHN1bShmbGFnMCkgPiAwCmlmKHN1bShmbGFnMCkgPiAwKXsKICBwcmludCgiV2FybmluZzogTm90IGFsbCBjb21iX2lkcyB3ZXJlIHF1YW50aWxlIG5vcm1hbGl6ZWQiKQp9ZWxzZXsKICBwcmludCgiQWxsIGNvbWJfaWRzIHF1YW50aWxlIG5vcm1hbGl6ZWQiKQp9CgphbnkoaXMubmEoZGF0YV9xbm9ybSkpCmBgYAoKIyMjIyBQbG90dGluZyBxbm9ybSBkYXRhIGZvciBRQyB2aXN1YWxpemF0aW9uCmBgYHtyfQojTWFrZSBib3hwbG90IG9mIGVhY2ggc2FtcGxlCmNvbF9tZXRhX3NlbGVjdCA9IGRwbHlyOjpzZWxlY3QoY29sX21ldGEsIGlkLCBjb21iX2lkKSAlPiUgZGlzdGluY3QoKQpkYXRhX3Fub3JtX2RmID0gYXMuZGF0YS5mcmFtZShkYXRhX3Fub3JtKQpkYXRhX3Fub3JtX2RmJGdlbmVfaWQgPSByb3duYW1lcyhkYXRhX3Fub3JtX2RmKQpkYXRhX3Fub3JtX2RmID0gdGlkeXI6OmdhdGhlcihkYXRhX3Fub3JtX2RmLCBrZXkgPSAiaWQiLCB2YWx1ZSA9ICJxbm9ybV92c3QiLCAtZ2VuZV9pZCkKZGF0YV9xbm9ybV9kZiA9IGRwbHlyOjpsZWZ0X2pvaW4oZGF0YV9xbm9ybV9kZiwgY29sX21ldGFfc2VsZWN0LCBieSA9ICJpZCIpCgpnZ3Bsb3QoZGF0YV9xbm9ybV9kZiwgYWVzKHggPSBpZCwgeSA9IHFub3JtX3ZzdCwgZ3JvdXAgPSBpZCkpICsKICBnZW9tX2JveHBsb3QoKSsKICBmYWNldF93cmFwKH5jb21iX2lkLCBzY2FsZXMgPSAiZnJlZV94IikrCiAgdGhlbWVfYncoYmFzZV9zaXplID0gMTQpCgpkYXRhX3Fub3JtX3BjYSA9IHJ1bl9wY2EoZGF0YV9xbm9ybSwgY29sX21ldGEgPSBjb2xfbWV0YSwgY29sb3JfY29sID0gInN0cmFpbl9pZCIsIHNjYWxlMCA9IFQpCgpgYGAKCmBgYHtyfQojU2F2ZSBxbm9ybSBkYXRhIHRvIGdjdApxbm9ybV9maWxlX3BhdGggPSBwYXN0ZTAob3V0ZGlyLCBwYXN0ZShydW5fbmFtZSxydW5fbmFtZV9zdWZmaXgsICJxbm9ybV9yZW1vdmVfbG93X2NvdW50X3NhbXBsZXNfYW5kX2dlbmVzIiwgc2VwID0gIl8iKSkKCnNhdmVfdG9fZ2N0KGRhdGFfbWF0ID0gZGF0YV9xbm9ybSwgZGF0YV9jZGVzYyA9IGNvbF9tZXRhLCBkYXRhX3JkZXNjID0gcm93X21ldGEsIHNhdmVmaWxlcGF0aCA9IHFub3JtX2ZpbGVfcGF0aCkKCmBgYAoKCgojIyBaLXNjb3JlIGFjcm9zcyBlYWNoIGdlbmUKIyMjIFByZXBhcmUgd2l0aCBRQwpgYGB7cn0KI0VhY2ggcm93IGlzIGEgZ2VuZSwgdmlzdWFsaXplIG1ldHJpY3MgdGhhdCB3aWxsIGJlIHVzZWQgZm9yIHotc2NvcmUgY2FsY3VsYXRpb24Kcm93X21lZGlhbnMgPSByb3dNZWRpYW5zKGRhdGFfcW5vcm0sIG5hLnJtID0gVCkKcGxvdCh4ID0gc2VxKDE6bGVuZ3RoKHJvd19tZWRpYW5zKSksIHkgPSByb3dfbWVkaWFucykKCnJvd19tYWRzID0gYXBwbHkoZGF0YV9xbm9ybSwgTUFSR0lOID0gMSwgRlVOID0gZnVuY3Rpb24oeCl7bWFkKHgsIG5hLnJtID0gVCl9KQpwbG90KHggPSBzZXEoMTpsZW5ndGgocm93X21hZHMpKSwgeSA9IHJvd19tYWRzKQoKYGBgCgojIyMgUGVyZm9ybSBaLXNjb3JlIGNhbGN1bGF0aW9uIHBlciByb3cKYGBge3J9CnpzY29yZV9saXN0ID0gbGlzdCgpCmRhdGFfenNjb3JlID0gZGF0YV9xbm9ybQoKZmxhZzAgPSByZXAoMCwgbGVuZ3RoKGNvbWJfaWRfdW5pcXVlKSkgI2ZsYWcgdG8gY2hlY2sgaWYgYWxsIGNvbWJfaWRzIHJhbgpmb3IobnVtIGluIDE6bGVuZ3RoKGNvbWJfaWRfdW5pcXVlKSl7CiAgICBpZHggPSB3aGljaChjb2xfbWV0YSRjb21iX2lkID09IGNvbWJfaWRfdW5pcXVlW251bV0pCiAgICBzdW0obGVuZ3RoKGlkeCkpCiAgICBpZihzdW0obGVuZ3RoKGlkeCkpPjApewogICAgICAgenNjb3JlX2xpc3RbW251bV1dID0gcm9idXN0X3pzY29yZShkYXRhX3Fub3JtWyxpZHhdLCBkaW0wID0gMSkgI3otc2NvcmUgYWxvbmcgcm93cywgb3V0cHV0cyB6c2NvcmVzLCBtZWRpYW5zIGFuZCBtYWRzCiAgICAgICBkYXRhX3pzY29yZVssaWR4XSA9IHpzY29yZV9saXN0W1tudW1dXVtbMV1dICNvbmx5IHRha2luZyB0aGUgei1zY29yZXMKICAgIH1lbHNlewogICAgICBmbGFnMFtudW1dID0gMQogICAgfQp9CgojQW55IGNvbWJfaWRzIGRpZG4ndCBnZXQgei1zY29yZWQ/CmlmKHN1bShmbGFnMCkgPiAwKXsKICBwcmludCgiV2FybmluZzogTm90IGFsbCBjb21iX2lkcyB3ZXJlIHN0YW5kYXJkaXplZCIpCn0KYGBgCgojIyMgQ2xpcCBleHRyZW1lIHotc2NvcmUgdmFsdWVzCmBgYHtyfQojQ2xpcCB6LXNjb3JlIHZhbHVlcyB0byBhY2NvdW50IGZvciBleHRyZW1lIGNhc2VzIHdoZXJlIGdlbmVzIHdpdGggbG93IGNvdW50IGhhZCBsb3cgdmFyaWF0aW9uIChpLmUuIGxvdyBtZWRpYW4gYWJzb2x1dGUgZGV2aWF0aW9uKQojTG9vayBhdCBnZW5lcycgei1zY29yZSBkaXN0cmlidXRpb24KaGlzdChkYXRhX3pzY29yZSwgYnJlYWtzID0gMTAwKQpxdWFudGlsZShkYXRhX3pzY29yZSwgYygwLjAwMSwgMC4wMDUsIDAuMDEsIDAuMDUsIDAuOTUsIDAuOTksIDAuOTk1LCAwLjk5OSkpCgpkYXRhX3pzY29yZSA9IGJhc2U6OnBtYXgoZGF0YV96c2NvcmUsIC0xKmNsaXBfdmFsdWUpCmRhdGFfenNjb3JlID0gYmFzZTo6cG1pbihkYXRhX3pzY29yZSwgY2xpcF92YWx1ZSkKCmBgYAoKCiMjIyBQbG90IFBDQQpgYGB7cn0KZGF0YV96c2NvcmVfcGNhID0gcnVuX3BjYShkYXRhX3pzY29yZSwgY29sX21ldGEgPSBjb2xfbWV0YSwgY29sb3JfY29sID0gInN0cmFpbl9pZCIsIHNjYWxlMCA9IFQpCmBgYAoKIyMjIFNhdmUgei1zY29yZXMgaW50byBHQ1QgZmlsZQpgYGB7cn0KaWYoIWFsbChjb2xuYW1lcyhkYXRhX3pzY29yZSkgPT0gY29sX21ldGEkaWQpKXsKICBwcmludCgiV2FybmluZzogY29sdW1uIG5hbWVzIG9mIG1hdHJpeCAmIG1ldGFkYXRhIElEcyBkbyBub3QgbWF0Y2giKQp9CmlmKCFhbGwocm93bmFtZXMoZGF0YV96c2NvcmUpID09IHJvd19tZXRhJGdlbmVfaWQpKXsKICBwcmludCgiV2FybmluZzogcm93IG5hbWVzIGRvIG5vdCBtYXRjaCBnZW5lIElEcyBpbiBtZXRhZGF0YSIpCn0KCnNhdmVfdG9fZ2N0KGRhdGFfbWF0ID0gZGF0YV96c2NvcmUsIGRhdGFfY2Rlc2MgPSBjb2xfbWV0YSwgZGF0YV9yZGVzYyA9IHJvd19tZXRhLCBzYXZlZmlsZXBhdGggPSBwYXN0ZTAob3V0ZGlyLCBwYXN0ZShydW5fbmFtZSxydW5fbmFtZV9zdWZmaXgsICJ6c2NvcmVfcmVtb3ZlX2xvd19jb3VudF9zYW1wbGVzX2FuZF9nZW5lcyIsIHNlcCA9ICJfIikpKQoKbnVtX3Jvd196c2NvcmUgPSBkaW0oZGF0YV96c2NvcmUpWzFdCm51bV9jb2xfenNjb3JlID0gZGltKGRhdGFfenNjb3JlKVsyXQoKYGBgCgojIyBDb2xsYXBzZSByZXBsaWNhdGVzCmBgYHtyfQp6c2NvcmUgPSBjbWFwUjo6cGFyc2VfZ2N0eChwYXN0ZTAob3V0ZGlyLCBwYXN0ZShydW5fbmFtZSxydW5fbmFtZV9zdWZmaXgsICJ6c2NvcmVfcmVtb3ZlX2xvd19jb3VudF9zYW1wbGVzX2FuZF9nZW5lcyIsIHNlcCA9ICJfIiksICJfbiIsIG51bV9jb2xfenNjb3JlLCAieCIsIG51bV9yb3dfenNjb3JlLCAiLmdjdCIpKQoKY29sX21ldGEgPSB6c2NvcmVAY2Rlc2MKZGF0YV96c2NvcmUgPSB6c2NvcmVAbWF0CgoKI2lkZW50aWZ5IHJlcGxpY2F0ZXMKI3BlcmZvcm0gbW9kZXJhdGVkIGNvbGxhcHNpbmcgb2YgYXZlcmFnZXMKY29uZGl0aW9uX2lkX3VuaXF1ZSA9IHVuaXF1ZShjb2xfbWV0YSRjb25kaXRpb25faWQpCm51bV9nZW5lcyA9IGRpbShkYXRhX3pzY29yZSlbMV0KbnVtX3NhbXBsZXMgPSBkaW0oZGF0YV96c2NvcmUpWzJdCm51bV9jb25kaXRpb25zID0gbGVuZ3RoKGNvbmRpdGlvbl9pZF91bmlxdWUpCgpkYXRhX3JlcF9jb2xsYXBzZSA9IG1hdHJpeChOQSwgbnJvdyA9IG51bV9nZW5lcywgbmNvbCA9IG51bV9jb25kaXRpb25zKQoKYWxsKGNvbF9tZXRhJGlkID09IGNvbG5hbWVzKGRhdGFfenNjb3JlKSkKCmZvcihudW0gaW4gMTpudW1fY29uZGl0aW9ucyl7CiAgaWR4ID0gd2hpY2goY29sX21ldGEkY29uZGl0aW9uX2lkID09IGNvbmRpdGlvbl9pZF91bmlxdWVbbnVtXSkKICBjaWR4ID0gd2hpY2goY29sbmFtZXMoZGF0YV96c2NvcmUpICVpbiUgY29sX21ldGEkaWRbaWR4XSkKICBkYXRhX3JlcF9jb2xsYXBzZVssbnVtXSA9IG1vZHpzX2NvbGxhcHNlX2F2ZyhkYXRhID0gZGF0YV96c2NvcmVbLGNpZHhdLCByaWR4ID0gc2VxKDEsIG51bV9nZW5lcywgYnkgPSAxKSkKfQpjb2xuYW1lcyhkYXRhX3JlcF9jb2xsYXBzZSkgPSBjb25kaXRpb25faWRfdW5pcXVlCnJvd25hbWVzKGRhdGFfcmVwX2NvbGxhcHNlKSA9IHJvd25hbWVzKGRhdGFfenNjb3JlKQoKYGBgCgoKCiMjIyBTYXZlIHJlcGxpY2F0ZSBjb2xsYXBzZWQgZGF0YSBhcyBnY3QgZmlsZQpgYGB7cn0KI1JldmFtcCBjb2x1bW4gbWV0YWRhdGEgc2luY2UgcmVwbGljYXRlcyBhcmUgY29sbGFwc2VkCmNvbF9tZXRhX2NvbGxhcHNlID0gc2VsZWN0KHNhbXBsZV9tZXRhLCBjb25kaXRpb25faWQsIHBlcnRfaWQsIHBlcnRfZG9zZSwgcGVydF9pZG9zZSwgcGVydF9kb3NlX3VuaXQsIHBlcnRfaW5hbWUsIHBlcnRfaXRpbWUsIHBlcnRfdGltZSwgcGVydF90aW1lX3VuaXQsIHBlcnRfdHlwZSwgcHJvamVjdF9pZCwgc3RyYWluX2lkLCB4X2Rvc2VfbWljX2hpZ2hfaW5vY3VsdW0pICU+JSBkaXN0aW5jdCgpCgpjb2xfbWV0YV9jb2xsYXBzZV9zdWIgPSBkYXRhLmZyYW1lKGNvbmRpdGlvbl9pZCA9IGNvbmRpdGlvbl9pZF91bmlxdWUpCmNvbF9tZXRhX2NvbGxhcHNlX3N1YiA9IGxlZnRfam9pbihjb2xfbWV0YV9jb2xsYXBzZV9zdWIsIGNvbF9tZXRhX2NvbGxhcHNlLCBieSA9ICJjb25kaXRpb25faWQiKQphbGwoY29sbmFtZXMoZGF0YV9yZXBfY29sbGFwc2UpID09IGNvbF9tZXRhX2NvbGxhcHNlX3N1YiRjb25kaXRpb25faWQpCmFsbChyb3duYW1lcyhkYXRhX3JlcF9jb2xsYXBzZSkgPT0gcm93X21ldGEkZ2VuZV9pZCkKCnNhdmVfdG9fZ2N0KGRhdGFfbWF0ID0gZGF0YV9yZXBfY29sbGFwc2UsIGRhdGFfY2Rlc2MgPSBjb2xfbWV0YV9jb2xsYXBzZV9zdWIsIGRhdGFfcmRlc2MgPSByb3dfbWV0YSwgc2F2ZWZpbGVwYXRoID0gcGFzdGUwKG91dGRpciwgcGFzdGUocnVuX25hbWUscnVuX25hbWVfc3VmZml4LCAibW9kenNjb3JlX3JlbW92ZV9sb3dfY291bnRfc2FtcGxlc19hbmRfZ2VuZXNfcmVwX2NvbGxhcHNlIiwgc2VwID0gIl8iKSkpCgpudW1fcm93X2NvbGxhcHNlID0gZGltKGRhdGFfcmVwX2NvbGxhcHNlKVsxXQpudW1fY29sX2NvbGxhcHNlID0gZGltKGRhdGFfcmVwX2NvbGxhcHNlKVsyXQpgYGAKCiMjIyBQbG90cyBhbmQgUUMKYGBge3J9ICAgIApjb2xfbWV0YV9jb2xsYXBzZV9zdWIkaWQgPSBjb2xfbWV0YV9jb2xsYXBzZV9zdWIkY29uZGl0aW9uX2lkCmRhdGFfenNjb3JlX3JlcF9jb2xsYXBzZV9wY2EgPSBydW5fcGNhKGRhdGFfcmVwX2NvbGxhcHNlLCBjb2xfbWV0YSA9IGNvbF9tZXRhX2NvbGxhcHNlX3N1YiwgY29sb3JfY29sID0gInBlcnRfZG9zZSIsIHNjYWxlMCA9IFQpCmBgYAoKCiMjIyBBZGQgaW5hY3RpdmUgaW5mbyBhbmQgZG9zZSByYW5rcyB0byB0aGUgY29sdW1uIG1ldGFkYXRhIGlmIGRlc2lyZWQgKG9wdGlvbmFsKQpgYGB7cn0KZGF0YV9yZXBfY29sbGFwc2VfcGF0aCA9IHBhc3RlMChvdXRkaXIsIHBhc3RlKHJ1bl9uYW1lLHJ1bl9uYW1lX3N1ZmZpeCwgIm1vZHpzY29yZV9yZW1vdmVfbG93X2NvdW50X3NhbXBsZXNfYW5kX2dlbmVzX3JlcF9jb2xsYXBzZSIsIHNlcCA9ICJfIiksICJfbiIsIG51bV9jb2xfY29sbGFwc2UsICJ4IiwgbnVtX3Jvd19jb2xsYXBzZSwgIi5nY3QiKQpkYXRhX3JlcF9jb2xsYXBzZSA9IHBhcnNlX2djdHgoZGF0YV9yZXBfY29sbGFwc2VfcGF0aCkKZGF0YV9yZXBfY29sbGFwc2VfbWF0ID0gZGF0YV9yZXBfY29sbGFwc2VAbWF0CgojRGV0ZXJtaW5lIHdlYWsgc2lnbmFsIHNlcGFyYXRlbHkgYmFzZWQgb24gREVTZXEyIHAtdmFsdWVzCmluYWN0aXZlcyA9IHJlYWRSRFMoJy9Vc2Vycy9zb3BoaWVfY2hlbi9kb3dubG9hZHMvcHNhX3JuYXNlcV9tYW51c2NyaXB0X3Jwcm9qZWN0LW1haW4vb3V0cHV0L21lb3dfdnN0X3JlbW92ZV9zYW1wbGVzX2FuZF9nZW5lcy5yZHMnKQoKI2NvbmRpdGlvbl9pZHMgY2hhbmdlZCB0byBpbmNsdWRlIDUgZGlnaXRzIGFmdGVyIHRoZSBkZWNpbWFsIGluIHRoZSBwZXJ0X2Rvc2UgdG8gYmUgY29uc2lzdGVudCB3aXRoIHByZXZpb3VzIHJ1bnMKaW5hY3RpdmVzJG51bWVyYXRvcl9jb25kaXRpb25faWRfb2xkID0gaW5hY3RpdmVzJG51bWVyYXRvcl9jb25kaXRpb25faWQKaW5hY3RpdmVzID0gbXV0YXRlKGluYWN0aXZlcywgbnVtZXJhdG9yX2NvbmRpdGlvbl9pZD1wYXN0ZShudW1lcmF0b3JfcHJvamVjdF9pZCwgbnVtZXJhdG9yX3BlcnRfaWQsIHBhc3RlMChzcHJpbnRmKCIlLjVmIiwgcm91bmQobnVtZXJhdG9yX3BlcnRfZG9zZSwgNSkpLCAidU0iKSwgbnVtZXJhdG9yX3N0cmFpbl9pZCwgIjkwbWluOkEiLCBzZXAgPSAiOiIpKQoKY2Rlc2MxID0gZGF0YV9yZXBfY29sbGFwc2VAY2Rlc2MKY29sbmFtZXMoY2Rlc2MxKQpjb2xuYW1lcyhpbmFjdGl2ZXMpCmFueSh0b2xvd2VyKGNkZXNjMSRjb25kaXRpb25faWQpICVpbiUgdG9sb3dlcihpbmFjdGl2ZXMkbnVtZXJhdG9yX2NvbmRpdGlvbl9pZCkpCmNkZXNjMSA9IG11dGF0ZShjZGVzYzEsIGluYWN0aXZlID0gaWZlbHNlKHRvbG93ZXIoY29uZGl0aW9uX2lkKSAlaW4lIHRvbG93ZXIoaW5hY3RpdmVzJG51bWVyYXRvcl9jb25kaXRpb25faWQpLCBUUlVFLCBGQUxTRSkpCmNkZXNjMSA9IGNkZXNjMSAlPiUKICBncm91cF9ieShzdHJhaW5faWQsIHBlcnRfaWQsIHBlcnRfdGltZSwgcGVydF9pdGltZSkgJT4lCiAgYXJyYW5nZShwZXJ0X2Rvc2UpICU+JQogIG11dGF0ZShkb3NlX3JhbmsgPSAxOm4oKSkgJT4lCiAgdW5ncm91cCgpCgphbnkoY2Rlc2MxJGlkICVpbiUgY29sbmFtZXMoZGF0YV9yZXBfY29sbGFwc2VAbWF0KSkKCiNyZW9yZGVyIGNkZXNjCmNpZHggPSBtYXRjaChjb2xuYW1lcyhkYXRhX3JlcF9jb2xsYXBzZUBtYXQpLCBjZGVzYzEkaWQpCmNkZXNjMSA9IGNkZXNjMVtjaWR4LCBdCmFsbChjZGVzYzEkaWQgPT0gY29sbmFtZXMoZGF0YV9yZXBfY29sbGFwc2VAbWF0KSkKc3RyKGNkZXNjMSkKCmRhdGFfcmVwX2NvbGxhcHNlQGNkZXNjID0gYXMuZGF0YS5mcmFtZShjZGVzYzEpCgp3cml0ZV9nY3QoZGF0YV9yZXBfY29sbGFwc2UsIHBhc3RlMChvdXRkaXIsIHBhc3RlKHJ1bl9uYW1lLCBydW5fbmFtZV9zdWZmaXgsICdtb2R6c2NvcmVfcmVtb3ZlX2xvd19jb3VudF9zYW1wbGVzX2FuZF9nZW5lc19yZXBfY29sbGFwc2VfaW5hY3RpdmVsYWJlbCcsIGRhdGUwLCBzZXAgPSAiXyIpKSkKCmBgYAo=